function f = UrbanEquilibriumFunction_poorneigh(pp,pop_C,pop_F,pop_P,xiH_F,xiH_P,xiH_C,XS,XM,varphi,QA_eq, alphaA,alphaM,alphaS,cAbar,cSbar, Xn,Xr,XrXr,XnXn,f_XrXn_C,f_XrXn_F,f_XrXn_P, L_r_M_only,L_n_M_only,Q_M_only,pSM_M_only);


   
%% unpackaging prices   
p_SU    = pp(1);
p_M     = pp(2);
   
%% constructing measures of skills distributions in Urban Regions (U= C + F)    
   f_XrXn_U   =  f_XrXn_C +f_XrXn_F+f_XrXn_P;
   
Xn_f_XrXn_C   =  XnXn.*f_XrXn_C;
Xn_f_XrXn_F   =  XnXn.*f_XrXn_F;
Xn_f_XrXn_P   =  XnXn.*f_XrXn_P;
Xn_f_XrXn_U   =  XnXn.*f_XrXn_U;

Xr_f_XrXn_C   =  XrXr.*f_XrXn_C;
Xr_f_XrXn_F   =  XrXr.*f_XrXn_F;
Xr_f_XrXn_P   =  XrXr.*f_XrXn_P;
Xr_f_XrXn_U   =  XrXr.*f_XrXn_U;

factor_wn  =  varphi^( varphi/(1-varphi) )*(1-varphi)*(XM/XS^varphi)^(1/(1-varphi));

if p_SU/p_M < pSM_M_only,
    wrU     =       varphi*p_M*XM*(L_n_M_only/L_r_M_only)^(1-varphi);
    wnU     =   (1-varphi)*p_M*XM*(L_n_M_only/L_r_M_only)^(-varphi);
else
    wrU                             =    p_SU*XS;
    wnU                             =    factor_wn*( p_M /(p_SU)^varphi )^(1/(1-varphi));  % these are the wages implied by the case when both S and M are produced at U
end

%% Individual Earnings
    yU                              =    max(wrU*XrXr,wnU*XnXn);
    
%% Consumption of Services
    ybar_F                          =    p_SU*cSbar*(1-alphaS)/alphaS + cAbar + p_M*xiH_F;        % income threshold for the consumption of Services, households in slums
    ybar_P                          =    p_SU*cSbar*(1-alphaS)/alphaS + cAbar + p_M*xiH_P;        % income threshold for the consumption of Services, households in poor neighborhoods
    ybar_C                          =    p_SU*cSbar*(1-alphaS)/alphaS + cAbar + p_M*xiH_C;        % income threshold for the consumption of Services, households in cities
    ICS_F                           =    zeros(size(XrXr));  % Indicator, positive consumption, households in slums
    ICS_P                           =    zeros(size(XrXr));  % Indicator, positive consumption, households in poor neighborhoods
    ICS_C                           =    zeros(size(XrXr));  % Indicator, positive consumption, households in cities
    ICS_F(yU>=ybar_F)               =    1;
    ICS_P(yU>=ybar_P)               =    1;    
    ICS_C(yU>=ybar_C)               =    1;
% Aggregation: Consumption of Services
    yTop_F                  =    trapz( Xr,trapz(Xn, yU.*f_XrXn_F.*ICS_F));    % Slums
    fTop_F                  =    trapz( Xr,trapz(Xn,     f_XrXn_F.*ICS_F));
    yTop_P                  =    trapz( Xr,trapz(Xn, yU.*f_XrXn_P.*ICS_P));    % Poor neighborhood
    fTop_P                  =    trapz( Xr,trapz(Xn,     f_XrXn_F.*ICS_P));
    yTop_C                  =    trapz( Xr,trapz(Xn, yU.*f_XrXn_C.*ICS_C));    % Cities
    fTop_C                  =    trapz( Xr,trapz(Xn,     f_XrXn_C.*ICS_C));

   
% Aggregation: Supplies of Skills and Employment    
    IOCC_r_U                        =    zeros(size(XrXr));  % indicator, occupation choices
    IOCC_r_U(wnU*XnXn < wrU*XrXr)   =    1;
    LrU                             =    trapz( Xr,trapz(Xn, Xr_f_XrXn_U .*IOCC_r_U));      %  effective labor in routine tasks, urban areas
    ErU=    trapz( Xr,trapz(Xn,    f_XrXn_U .*IOCC_r_U));      %  employment in routine tasks, urban areas
    LnU=    trapz( Xr,trapz(Xn, Xn_f_XrXn_U .*(1-IOCC_r_U)));  %  effective labor in non-routine tasks, urban areas
    EnU=    trapz( Xr,trapz(Xn,    f_XrXn_U .*(1-IOCC_r_U)));  %  employment in non-routine tasks, urban areas
if p_SU/p_M<pSM_M_only,
    QM=   Q_M_only;
    QSU=   0; 
else
    LrM                                       =   (varphi*XM*p_M/XS/p_SU)^(1/(1-varphi))*LnU;
    QM  =   XM*(LrM)^varphi*(LnU)^(1-varphi);
    QSU =   XS*(LrU-LrM);          
end
 

    
    
%% Equilibrium Pricing Conditions
    PP_SU1  =    alphaS*( yTop_F+yTop_P+yTop_C - fTop_C*(cAbar+p_M*xiH_C) - fTop_F*( cAbar+p_M*xiH_F ) - fTop_P*( cAbar+p_M*xiH_P ) )/( QSU +(1-alphaS)*cSbar*(fTop_F+fTop_P+fTop_C) );
    PP_MM1  =    alphaM/alphaA*( QA_eq - cAbar )/( QM -xiH_C*pop_C-xiH_F*pop_F-xiH_P*pop_P );
    
 
    
%% Distance (Euclidean)
%f  =  (i_fr-fr)^2+(i_m_y_r-m_y_r)^2 +(i_m_y_n-m_y_n)^2 +(i_s_y_r-s_y_r)^2 +(i_s_y_n-s_y_n)^2;

 f = (PP_SU1-p_SU)^4 + (PP_MM1-p_M)^4;  
 
 
